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Abstract. A method is described for embedding a deformable, elastic, membrane 
within a lattice Boltzmann fluid. The membrane is represented by a set of massless 
points which advect with the fluid and which impose forces on the fluid which are 
derived from a free energy functional with a value which is dependent upon the 
geometric properties of the membrane. The method is validated in two dimensions 
with a free energy functional which imposes the constraint of constant membrane 
length, constant enclosed bending rigidity and a preferred curvature. The 

method is shown to recover the expected equilibrium shape in the absence of flow 
and deformation in the presence of an applied shear flow. The method may have 
applications in a number of mesoscopic simulations, including discrete models of blood 
cells. 
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1. Introduction 

A number of biological objects can be represented as vesicles formed by polymerised 
membranes (Alberts et al 2002, Lipowsky 1991). While the conformations of such 
vesicles are determined by the elastic properties of the membrane (Lipowsky 1991), 
their dynamics in a flowing fluid is altered by the flow of the fluid outside and inside 
vesicles. Both factors should be taken into account in simulation of the membranes 
immersed into fluid host. 

The lattice Boltzmann (LB) method (Succi 2001) provides a convenient method for 
introducing fluid flow in the presence of boundaries. An extension of the LB method 
to model fluid membranes in which the molecules rapidly diffuse within the membrane 
has been reported (Stelitano and Rothman 2000). However, this method cannot be 
straightforwardly generalised to the polymerised membranes since it does not allow for 
the extensional elastic properties which would be associated with such a membrane. An 
alternative approach is to model the membrane as a geometric object immersed into the 
LB fluid and this is the method adopted in this paper. A similar approach has been 
developed for the simulation of a polymer chain in an LB solvent (Ahlrichs and Diinweg 
1998). 

The aim of the present paper is to develop an LB method for the polymerised 
membranes. The main purpose of the paper is methodological and we therefore 
demonstrate the effectiveness of the method for a simple two-dimensional case. However, 
we discuss the generalisation to three dimensions in section 0D It is important to note 
that the motivation for the work in this paper is ultimately to develop models for flows 
in which there are embedded a large number of deformable membranes; this will form 
the basis of modelling blood flow at the veinule scale. For this reason, the finer details 
of the membrane properties are ignored (eg membrane dissipation). 

The paper is organised as follows. In section 121 the method for introducing the 
membrane into the LB scheme is described. The explicit expressions for the forces 
arising from the membrane are derived in section El The results of the simulations of 
the relaxation of a closed membrane to its equilibrium state without and with shear 
flow are presented in section HJ Conclusions and possible extensions and applications 
are discussed in sectional 

2. Lattice Boltzmann 

The basics of the LB method have been described in the literature (Succi 2001). The 
fluid in the LB method is considered as a field of the population densities /j(r, t), which 
indicate the amount of fluid present at the lattice site r at the discrete time t and 
moving with the velocity q associated with the i-th lattice direction. The models with 
n velocities on a simple cubic lattice of dimension d are usually referred to as DdQn. 
The LGBK algorithm may be represented by the equation 

fi(r + *6 t , t + S t ) = fa, t) + -(fP - fi) + G i} (1) 
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where S t represents the time step, r controls the kinematic shear viscosity of the lattice 
fluid through the relation 
2r - 1 



v 



6 

and the 'forcing' term 



St, (2) 
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may be used to impress an external force F M , where c s is the speed of sound and the U 
are determined to achieve isotropy of the fourth-order tensor of velocities and Galilean 
invariance. Note that the expressions (|Tfl-([3*J) are independent of the spatial dimension. 
Velocity moments give the lattice fluid's density and momenta through 

P = E/* = Eif, (4) 

i i 

pv = Y,c l f l = Y,^°\ (5) 



where the equilibrium distribution function is 



fi (j>,v) = Up 



v • c, |v| 2 (V • c, 
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The form of the equilibrium distribution function (JBJ) ensures that relations ()4I5|) are 
satisfied and also determines the nonviscous part of the momentum-flux tensor of the 
lattice fluid, 

n i°/3 = E fx C iaCi/3 = C 2 s pS a(3 + fWaVp. (7) 
i 

The membrane is represented by a discrete set of points corresponding to the equidistant 
values of a parameter s. The difference As = Si+i — Sj between the values of the 
parameter s for the consecutive points is chosen to be comparable with the distance 
between LB nodes. 

The membrane is treated purely as a geometric object and, contrary to (Ahlrichs 
and Diinweg 1998), the points of the membrane always move with the velocity of the 
underlying LB fluid determined from equation (J3J). Generally, the position of the points 
does not coincide with the location of LB nodes, so a weighted average is used based on 
the velocities at the nodes which bound the lattice primitive cell in which the point is 
situated. The weight of the contribution from each node is taken to be proportional to 
the distance of the node from the membrane point under consideration. The geometric 
properties of the membrane determine the force that is in turn applied to the LB fluid 
according to equation (jBJ). Again, as the location of the point generally does not coincide 
with nodes, the force is distributed amongst the same set of nodes in accordance with 
their distance from the point at which the force originates. The value of the force is 
determined in section El 

It should be noted that a membrane represented in this way will not strictly 
conserve its internal mass, and will therefore be slightly permeable. However, since 
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we are modelling polymerised membranes which are themselves weakly permeable, this 
is not considered to be a serious defect in the approach. 

3. Forces 

We describe the shape of the membrane by the vector function r(s) having the 
components x(s) and y(s). We choose the parameter s so that for the undeformed 
membrane (without internal stress) s coincides with the arc length parameter / of the 
curve. The parameter s spans the interval (0,L ), L being the equilibrium length of 
the membrane. In general, the length of the membrane is determined by the expression 

L = [ dl = [ L °u(s)ds. (8) 



* L JO 

Here dl = u(s)ds is the length element, 

u(s) = ^x'\s)+y%s). (9) 

If the membrane is stretched/compressed, the free energy increases. The excess free 
energy is 

A L [v{s)\ = % l L ° («(*) - if ds, (10) 
2 Jo 

a being the membrane compressibility. The free energy arising from the bending 
elasticity of the membrane can be taken in form (Canham 1970, Helfrich 1973, Evans 
1974): 

A K = ^J L (K(r)-K fdl, (11) 

K(r) being the curvature, K being the preferred curvature and k being the bending 
rigidity coefficient. The curvature can be represented as a function of the parameter s 

x'( S )y"(s)-y'(s)x"( S ) 

K{S) = ^ • (12) 

and hence we can write 

M*(*)\ = ? f L ° ( K ( s ) - K o? »(8)d8. (13) 
I Jo 

The contribution to the free energy due to the surface tension is A s = aL, or 

A s = a [ ° u(s)ds, (14) 
Jo 

a being the surface tension coefficient. The fluid is assumed to be compressible and the 
equilibrium two-dimensional 'volume' of the fluid enclosed by the membrane is taken to 
be Vq. The excess pressure is 

p=-P(V[r(s)]/V -l), (15) 

where (3 is the fluid compressibility, and the 'volume' of the droplet V[r(i)] is the 
functional of the membrane shape 

V[r(s)} = [ L ° y(s)x'(s)ds. (16) 
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JL(V[r(s)]-V ) 



The free energy due to the fluid compressibility is 

A v = ( V pdV (17) 

JVo 

or, after integration, 
A v [r(s)] 

As a result, the free energy of the interface can be represented as a following 
functional of the membrane shape: 

A[r(s)} = A L [r(s)} + A K [r(s)] + A s [r(s)\ + A v [v(s)\, (19) 

where y4 L [r(s)], Ax[r(s)], A^-[r(s)] and Ay[r(s)] are given by formulas (fTU|) . (|TH|). (JT3J), 
and (fTHj) . correspondingly 

The force F from the element dl of the membrane is found as the variational 
derivative of the free energy, F(s) = — 5A[r(s)]/5r(s), and can be represented in form 



F(s) = F L (s) + F K (s) + F s (s) + F v (s) 
with the components 

F L (s) = -a (K(s)n(s) + d 2 r/ds 

1 



Fk(s) 

F V (S) 



K 
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d 2 K(s) 
dl 2 



n(s), 



-(jii'(s)n(s), 
-/3(V[r( a )]/V 



l)n(a) 



(20) 

(21) 

(22) 

(23) 
(24) 



corresponding to the membrane compressibility, bending elasticity, surface tension and 
fluid compressibility. In formulas ()21fj2"l"j) . n(s) is the unit vector normal to the 
membrane with the components n x (s) = y'(s)/u(s) and n y (s) = —x'(s)/u(s). The 
expression for the force (|22j) is equivalent to that in (Stelitano and Rothman 2000) 
after it has been noted that there is an error in the quoted result which arises because 
the authors do not correctly account for the change in the metric tensor of the surface 
during the minimisation of the free energy (Lishchuk and Care 2005). Equation ()22j) 
also includes an additional contribution due to the preferred curvature K . 

The functions x(s) and y(s) can be approximated by polynomials. We use the 
second order polynomials, 

x(s) = a + aiS + a 2 s 2 , 
y(s) =b + b v s + b 2 s 2 . 

In order to determine the coefficients and b^ {k = 0,1,2} at the z-th point of 
the discretised membrane, we require the values of the functions at this and two 
neighbouring points to coincide with the corresponding positions, 

1 ) %i— 1 

x(si) = Xi , (25) 

^(■Si+l) = 
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This represents a system of equations for a^, and there is an analogous system for 
The values of these coefficients give the explicit form of the functions x(s) and y(s) in 
the vicinity of each point s, that can be used to find its derivatives up to the second 
order. We note that the force on a point is different from the force on an element dl by 
a factor u(s). 

4. Validation 

The method described in this paper is intended as methodological; there are no 
experimental results available for a two dimensional system. The results we give below, 
demonstrate that the method behaves in a manner which is consistent with the expected 
behaviour of a two dimensional cell. The possible extension of the method to three 
dimensions is discussed in section El 

The simulations in this section were run on D2Q9 100 x 100 lattice with the periodic 
boundary conditions. The basis velocity vectors of the D2Q9 lattice and corresponding 
values of U are presented in the table ^ The primitive cell used for averaging the 
velocities, and distributing the forces, was taken to be a primitive square cell of the 
lattice whose corners are lattice nodes. The value of the relaxation time r = 0.8 was 
used in the simulations. The surface tension coefficient and the preferred curvature of 
the membrane were set to zero. The equilibrium distance between the membrane points 
As = 3.2 was used, and a typical membrane included 50 — 100 points. 

Table 1. The basis velocity vectors of the D2Q9 lattice and corresponding values of 



u. 
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a 


u 





(0,0) 


4/9 


1 


(1,0) 


1/9 


2 


(0,1) 


1/9 


3 


(-1,0) 


1/9 


4 


(o,-i) 


1/9 


5 


(1,1) 


1/36 


G 


(-1,1) 


1/36 


7 


(-1,-1) 


1/36 


8 


(1,-1) 


1/36 



Figure ^ shows the equilibrium shapes of the vesicle for different values of the 
parameter Q defined as 

Q = L /L c , (26) 

L c being the length of the circular membrane with the same enclosed area. The 
parameters are: k = 0.05, a = 0.01, (3 = 0.2. To provide the slight initial asymmetry, 
the initial shape of the droplet was an ellipse with half-axes 20 ± 0.2. 
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It should be noted that the same equilibrium shapes can be obtained without the 
including the effect of the embedding fluid; this was confirmed as one of the validations 
of the simulation. However, the inclusion of the LB fluid is necessary in order to recover 
the dynamics of the relaxation to the equilibrium shape, as is shown in Figure 121 for 
Q — 1.4 and different values of the bending rigidity k. The value of k also influences 
the rate of relaxation which is demonstrated in Figure El by the time dependence of the 
mean square velocity of the LB fluid for different values of k. 

To investigate the behaviour of the membrane in the flow, a simulation was 
undertaken in which shear flow was applied to an initially spherical membrane. Figure^] 
depicts the time evolution of the shape of the droplet, and Figure El shows the velocity 
field of the LB fluid. Apart from the imposed shear, the parameters of the simulation 
are the same as in Figure H with the parameter Q was taken to be equal to 1.4. In 
the final steady state the membrane rotates with the fluid; however, this behaviour is 
expected only to occur in a two-dimensional system. 

5. Conclusion 

A method has been described for embedding a deformable membrane into a LB fluid 
and results presented which validate the approach in two dimensions. The method does 
not take into account thermal fluctuations. If the fluctuations are small they simply 
result in a renormalisation of the bending rigidity coefficient (Palmer and Morse 1996), 
and no further modification of the method is necessary. Large fluctuation can be taken 
into account by adding the random stress to the LB algorithm (Cates et al 2004). 

The method could be generalised to three dimensions by employing the expression 
for the force due to bending rigidity of the two-dimensional membrane based on a 
corrected version of the result derived by (Stelitano and Rothman 2000) (see comment 
after equation (}2"4"|) ). and the method for calculating the curvature of the triangulated 
surface which has been developed by (Hamman 1993). Note that a grid would need to be 
created for the equilibrium shape of the membrane prior to LB simulation. This could 
be achieved by the direct numerical minimisation of the free energy of the membrane. 
After the minimisation, the equilibrium distances between points would be changed and 
the details of the area contributions would therefore need to be modified appropriately; 
work is currently in progress to implement such a three dimensional scheme. 

One possible application of the technique described in this paper is a more accurate 
representation of the flow of blood cells in confined geometries and in which membrane 
elasticity effects are taken into account. 

References 

Ahlrichs P and Dunweg B 1998 Int. J. Mod. Phys. C 9 1429-38 

Alberts B et al 2002 Molecular Biology of the Cell (New York: Garland Science) p 1549 

Canham P B 1970 J. Theor. Biol. 26 61-81 

Cates ME et al 2004 J. Phys.: Condens. Matter 16 S3903-15 



A deformable elastic membrane embedded in a lattice Boltzmann fluid 



8 



Evans E 1974 Biophys. J. 14 923-31 
Hamman B 1993 Computing Suppl. 8 139-53 
Helfrich W 1973 Z. Naturforsch. 28c 693-703 
Lipowsky R 1991 Nature 349 475-81 

Lishchuk S V and Care C M 2005 Phys. Rev. E submitted 
Palmer K M and Morse D C 1996 J. Chem. Phys. 105 11147-74 
Stelitano D and Rothman D H 2000 Phys. Rev. E 62 667-80 

Succi S 2001 The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (New York: Oxford 
University Press) p 288 



A deformable elastic membrane embedded in a lattice Boltzmann fluid 



9 




Q=1.4 Q=1.6 



Figure 1. Equilibrium shapes for different values of the parameter Q defined by J2HJ)- 




k=0.015 k=0.050 



Figure 2. Dynamics for different values of the bending rigidity k. 
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Figure 3. Time dependence of the mean square LB velocity for k = 0.025 (curve 1), 
k = 0.05 (curve 2), n = 0.1 (curve 3). 




Figure 4. Dynamics of the membrane in shear flow. 
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Figure 5. The velocity field of the LB fluid under shear. 



